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SIMULATION  OF  SOME  AUTOREGRESSIVE  MARKOVIAN  SEQUENCES 


OF  POSITIVE  RANDOM  VARIABLES 

by 

A.  J.  Lawrance 
and 

P.  A.  W.  Lewis 

ABSTRACT 

Methods  for  simulating  dependent  sequences  of  continuous 
positive-valued  random  variables  with  exponential,  Gamma,  mixed 
exponential  and  uniform  marginal  distributions  are  given.  In 
roost  cases  the  sequences  are  first-order,  linear  autoregressive, 
Markovian  processes.  A  two-parameter  family  of  this  type 
with  exponential  marginals  is  defined  and  its  transformation 
to  a  similar  multiplicative  process  with  uniform  marginals 
is  given.  It  is  shown  that  for  a  subclass  of  this  two- 
parameter  family  extension  to  mixed  exponential  marginals  is 
possible,  giving  a  model  of  broad  applicability  for  analyzing 
data  and  modelling  stochastic  systems.  Efficient  simulation 
of  some  of  these  schemes  is  discussed. 


1 .  INTRODUCTION 


In  a  recent  series  of  papers  [1,2,3,4,5,6,7,81  some 
simple  models  have  been  derived  for  stationary  dependent 
sequences  of  positive,  continuous  random  variables  with  given 
first-order  marginal  distributions.  In  general  the  dependency 
structure,  as  measured  by  second-order  joint  moments  (serial 
correlations)  mimics  that  of  the  usual  linear  mixed  auto- 
regressive-moving  average  (ARMA)  models  which  have  been  used 
for  so  long  in  time-series  analysis.  In  the  ARMA  models, 
which  are  defined  quite  generally,  there  is  in  usage  an 
implicit  assumption  of  marginal  normality  of  the  random  vari¬ 
ables.  This  is  clearly  not  the  case  if  the  random  variables 
are  positive,  say  the  times  between  events  in  a  series  of 
events  [9]  or  the  successive  response  times  at  a  computer 
terminal.  Thus  the  new  models  are  derived  to  accommodate 
situations  in  which  the  dependent  random  variables  have,  for 
instance,  exponential,  Gamma,  Laplace  and  mixed  exponential 
marginal  distributions.  The  exponential  case  is  the  most 
highly  developed,  with  the  nomenclature  [4]  EARMA (p,q) 
(exponential  process  with  mixed  moving  average-autoregressive 
structures  of  orders  p  and  q  respectively)  and  NEARMA(p,q) 

(new  EARMA (p,q) ) . 
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The  development  of  the  probabilistic  properties  of 
these  processes  is  given  in  the  referenced  papers,  applications 
to  queueing  models  and  computer  system  modelling  in  [10,11] 
while  development  of  estimation  and  testing  procedures  has 
just  begun. 

The  object  of  the  present  paper  is  to  define  and 
discuss  the  simulation  of  the  processes  on  digital  computers, 
though  for  the  sake  of  brevity  only  the  first-order  Markovian, 
autoregressive  case  is  considered.  The  simplicity  of  structure 
of  these  models — in  general  they  are  linear  additive  mixtures 
of  random  variables— makes  them  ideal  for  this  purpose .  How¬ 
ever  stationarity  conditions  are  sometimes  difficult  to  derive 
analytically  and  in  some  cases  it  is  not  simple  to  generate 
the  innovation  random  variables  in  the  processes.  A  striking 
example  of  this  is  the  case  of  the  Gamma  first-order  autore¬ 
gressive  process  for  which  an  efficient  means  of  simulation 
has  only  recently  been  found  [7]  for  some  parametric  values. 

It  is  shown  that  a  simple  transformation  of  the  exponential 
sequences  gives  a  direct  multiplicative  method  for  generating 
dependent  processes  with  uniform  marginals.  These  could  be 
the  basis  in  simulations  for  many  other  types  of  dependent 
sequences . 
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2.  EXPONENTIAL  AUTOREGRESSIVE  MARKOVIAN  SEQUENCES 

We  give  here  three  methods  of  generating  first-order 

autoregressive,  Markovian  sequences  with  exponential  marginal 

distributions.  The  first  two  are  defective  in  terms  of  their 

sample  path  properties  (the  first  more  so  than  the  second) 

while  the  third,  NEAR(l),  is  satisfactory  in  this  respect. 

The  defect  of  the  first  two  models  is  also  highlighted  by 

the  simulation  procedures  used;  they  can  be  generated  from 

one  sequence  of  exponential  variables. 

Note  that  autoregression  in  the  context  of  a  stochastic 

sequence  (Xn)  is  vaguely  used.  In  the  first  place  linear, 

additive  autoregression  is  usually  implied.  In  the  second 

place  autoregression  can  mean  that  in  the  defining  equation 

for  Xn  the  previous  value  enters  explicitly,  but  more 

particularly  it  means  that  the  conditional  expectation  of 

X_,  given  X_  .  =  x^  . ,  is  an  additive  linear  function  of 
n  n-±  n-± 

"h-i'* 


E(Xn|xn-l 


=  x^  . )  =  a  +  bx 
n-i 


n-1  * 


(1) 


The  Markovian  property  (first-order)  means  that  the  probability 
structure  of  Xn,  Xn+1,...  ,  given  Xn-1  =  xR  is  independent 

of  xn-2'  Xn-3'***  * 
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2A.  The  Exponential  DAR(l)  Process 

A  very  simple  exponential  autoregressive  Markovian 
sequence  is  generated  by  the  equation 


xn  -  V»-l  +  (1  - 


v_)e_ 

n  n 


(2) 


where  p{vn*l)  »  1  -  p{vn=0)  =  p  and  En,  n  =  1,2,...  are 
independent  exponential  random  variables  with  parameter  X; 


P{En  _<  x)  =  1  -  e~*x,  x  ^  0,  X  >  0 

=  0  ,  x  <  0 . 


(3) 


For  this  process  the  serial  correlations 
are 


pk  »  oorr(Xn,xn+k) 


and 


(4) 


E‘Wl  -  Vl>  =  plxn-l  +  (1  -  pl>/x-  <5» 


This  process  is  an  exponential  version  of  the  DAR(l)  process 
[12,13]  but  is  rather  useless  for  modelling  real  data  because 
simulations  of  the  process  show  runs  of  xn's  with  the  same 
value.  These  occur  when  xn_^  is  picked  successively  in  (2) , 
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rather  than  the  innovation  En.  Moreover  the  lengths  of  the 
runs  of  similar  values  are  geometrically  distributed. 

2B.  The  Exponential  EAR(l)  Process 

Another  model  is  derived  from  the  usual  linear  model 


Xn  *  pXn-l  +  en 


(6) 


in  which  the  i.i.d.  innovation  process  {en>  is  chosen  so 
that  the  xn'8  are  marginally  exponential (X) .  Gaver  and 
Lewis  [1]  show  that  for  this  to  be  true,  one  must  have 
0  <  p  <  1  and 


en  =  E 
n  n 


-  0 


w.p.  1-p, 
w.p.  p  , 


(7) 


where  (En),  as  throughout  the  paper,  are  i.i.d.  exponential (X) 

Again  -  p  and  E(xnlxn_i  ■  xn_i^  =  pixn-l  + 
as  at  (4)  and  (5)  for  the  exponential  DAR(l)  model.  The 
difference  is  in  the  sample  paths;  for  the  EAR ( 1 )  process 
simulations  show  runs  of  xn's  decreasing  geometrically. 

These  occur  when  only  Pxn_^  is  picked  in  (6) .  Again  the 
lengths  of  these  runs  are  geometrically  distributed. 


The  Markov  property  of  the  two  sequences  implies 

that  if  Xq  is  chosen  to  be  Eq,  an  exponential (X)  random 

variable  independent  of  E^  E2,  ...  ,  then  ,  X2,... 

forms  a  stationary  sequence. 

Naive  inspection  of  the  defining  equations  (2) ,  (6) 

and  (7)  suggest  that  to  generate  a  stationary  sequence  of 

length  N,  X^,...,  Xjj,  (N+l)  i.i.d.  exponential  deviates 

and  N  uniform  variates  (for  the  selection  process)  are 

needed.  However,  the  sequences  can  be  generated  from  only 

one  exponential  sequence;  this  is  possibly  related  to  the 

degeneracy  in  the  processes.  This  method  uses  the  memory less 

property  of  exponential (X)  variables,  namely  that  if  Er 

is  given  to  be  greater  than  a  constant  y,  then  En  "  Y 

is  again  exponential (X) . 

Thus  the  algorithm  is  to  initialize  by  setting 

XQ  =  Eq;  subsequently  set  XR  =  PXn-1  if  En£xp  = -£n(l-p)/X; 

otherwise  set  X  =  pX  .  +  (E  -  x  ) .  This  uses  the  fact 

n  n-l  n  p 

that,  from  (3),  P(En  £  xp)  =  p. 

Even  greater  efficiency  can  be  obtained,  though  this 

must  be  qualified  by  considerations  as  to  whether  the  X  's 

n 

are  to  be  generated  one  at  a  time  or  in  an  array;  whether 
a  subroutine  is  available  to  generate  exponential  random 
variables  faster  than  can  be  done  by  taking  logarithms  of 
uniform  deviates,  and  the  relative  speed  of  division  and 
generation  of  uniform  deviates. 
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The  more  efficient  scheme  recycles  uniform  variables. 


i.e.  if  0  is  given  to  be  between  constants  a  and  b, 
where  0  <  a  <  b  £  1,  then  (U-a)/(b-a)  is  a  uniform 
random  variable.  (Note  that  its  value  is  not  given,  only 
that  it  is  in  (a,b).)  Thus  to  generate  an  array 
X^,...,  *n'  of  dependent  exponential  variates  with  mean  1 
from  the  EAR(l)  process  we  have 

Algorithm  1  (EAR(l)  process) 

1.  Generate  U  and  set  Xq  -An  U 

2.  Generate  U  and  set  Y  *■  U 

3.  n  ♦  1 

4.  If  Y  £  p  go  to  7 

5.  ♦  PX^  .  -  An (Y-p)  +  An(l-p) 

n  n-i 

6.  Generate  U,  set  Y  U  and  go  to  8 

7.  X„  •*-  pX  .  ;  Y  Y/p 

n  n-i 

8 .  Set  n  *■  n+1 

9.  If  n  <_  N  go  to  4 

10.  Otherwise  exit. 

The  expected  number  of  uniform  deviates  required  in  this 
algorithm  is  1  +  (l-p)N,  which  is  less  than  the  number  N 
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2C.  The  Exponential  NEAR  (1)  Process 

A  broader  two-par ame ter  exponential  sequence  which 
is  a  first-order  autoregressive,  Markovian  process  and  an 
additive  linear  mixture  of  random  variables  is  given  by 
Lawrance  [7]  and  developed  by  Lawrance  and  Lewis  [5] .  Called 
NEAR(l) ,  the  sequence  is  defined  as 


pxn_!  W*P*  a 


X  =  e  + 
n  n 


n  =  1, 2, . . .  , 


w.p.  1- a 


(8) 


where  0  _<  a  <  1  and  0  £  8  <,  1.  It  can  be  shown  that  for 
the  XR  to  be  marginally  exponential (X )  the  innovation 
variable  en  must  be  generated  from  an  En  by  the  exponential 
mixture 


e 


n 


(1-a)  0En 


w.p. 


1~8 

1- (1-a) 8 


w.p. 


c<8 

l-fl-a)8 


n  1,2,...  (8) 


providing  a  and  $  are  not  both  equal  to  one.  When  a  =  0 
or  M  0  the  {XR}  are  exponential  i.i.d.,  whereas  when 
a  *  1  the  EAR(l)  model  given  at  (6)  and  (7)  is  obtained.  In 
fact  fixing  either  a  or  8  and  varying  the  other  parameter 
gives  an  exponential  model  with  a  full  positive  range  of 
serial  correlation  of  order  one,  since  it  is  easily  shown 
that 
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Again 


E  (X  |X  ,  =  x  .)  =  a0x  .  +  (l-a0 )  /X 
n  n-x  n-x  n~x 

=  Plxn-l  +  U-Pj.)/*  (ID 

and  Xq  =  Eq  gives  a  stationary  sequence.  The  NEAR(l)  pro¬ 
cess  allows  one  to  model  a  broader  class  of  exponential 
sequence  as  measured  either  by  sample  path  behavior  or  higher- 
order  joint  moments;  see  [5]  for  details. 

A  particularly  simple  case  occurs  when  0  *  1;  this 
model,  called  TEAR(l) ,  is  very  tractable  analytically  and, 
as  will  be  shown  below,  extends  easily  to  the  case  of  mixed 
exponential  distributions  for  the  XR. 

Note  that  in  the  NEAR(l)  process  the  innovation  e 

n 

is  always  present  unless  a  =  1  and  it  is  therefore  not 
possible  to  simulate  the  stationary  process  with  less  than 
N+l  uniform  variates.  The  simplest  method  to  generate 
dependent  exponential  variates  with  mean  1  from  the  NEAR(l) 
process  seems  to  be  the  following: 

Algorithm  2  (NEAR(l)  process) 

1.  Generate  U;  set  Xq  «-  -in  U;  y  «-  1-a; 

5  «-  (l-0)/[l  -  (l-a)0  ] . 

2.  Generate  U;  set  n  *■  1. 


4.  Otherwise  Y  ♦  (U  -  Y)/(l  -  Y) . 


5 .  If  Y  <  6  set  Xn  «-  3  Xn_1  -  An  Y  +  in  6  and  go  to  9 . 

6.  Otherwise  set  xn  «-  6xn_^  "  y6  An  [  (Y-fi  )/(l-6)  ]  and 
go  to  9. 

7.  If  Y  £  6  set  .  XR  -An  Y  +  An  6  and  go  to  9. 

8.  Otherwise  set  Xn  «•  An  [  (Y-6  )  /(1-6 )  ]  . 

9.  Set  n  n+1. 

10.  If  n  <_  N  generate  U  and  go  to  3.  Otherwise  exit. 

Note  that  for  stationary  array  of  N  Xn's,  exactly 
N+1  uniforms  are  required  and  therefore  it  could  be  advantageous 
to  generate  these  in  an  array  which  would  be  replaced  one  at 
a  time  by  the  Xn's.  Care  must  be  taken  with  the  recycling 
of  the  uniform  variates  U  if  y  =  l-*a  is  close  to  one  or 
zero.  In  that  case  it  is  probhbly  better  for  computational 
reasons  to  use  2 (N+1)  uniform  variates.  Note  that  y  *  1 
gives  the  EAR(l)  process. 


3.  UNIFORM  MARKOVIAN  SEQUENCES 

It  is  convenient  to  have  dependent  sequences  of  random 
variables  with  marginal  distributions  other  than  exponential. 
Before  discussing  other  solutions  to  the  Equation  (8)  we 
show  that  a  simple  transformation  of  the  NEAR(l)  process  gives 
a  two-parameter  family  of  Markovian  random  variables  with 


uniform  marginal  distributions.  It  is  well-known  that  an 
exponential  transformation  of  a  unit  exponential  random  vari¬ 
able  gives  a  uniformly  distributed  random  variable.  Thus  we 
have  from  (8)  and  (9)  the  multiplicative  model  for  a  uniform 
Markovian  sequence  {Xn>,  n  =  1,2,...  ; 


n 


e  • 
n 


w.p.  a 


=  e 


n 


w.p.  (1-a) 


n  1,2,... 


where 


(12) 


n 


*  u 


n 


_  j.(l-a)  8 

un 


w-p-  i^nrlrg 

a3 

w*p*  r-(l-a) 3 

n  =  1,2,... 


(13) 

(14) 


for  Un,  n  =  1,2,...,  i.i.d.  uniformly  distributed,  providing 
that  a  and  8  are  not  both  equal  to  one.  Again  if  Xg 
is  uniformly  distributed  and  independent  of  U^,  U ...  the 
sequence  is  stationary. 

An  algorithm  for  generating  this  two-parameter  uniform 

sequence  is  easily  adapted  from  Algorithm  2.  It  remains  to 

find  the  correlation  structure  and  the  regression  of  Xn 

on  X.  .. 
n-i 

To  do  the  former,  let  X^  be  a  NEAR(l)  sequence  with 
X  *  1,  so  that  the  sequence  XR  at  (12)  is  given  by 
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^  ■  e*p{-X£} .  Now  the  joint  Laplace-Stieltjes  transform  of 
X*,  xn-k  by  Lawrance  and  Lewis  [5]  as 


(s,  t) 


»  E{exp[-sX*  -  tx*  .]} 

n  n-k  (is) 

v  k  .  k  k-1  .  i 

-  o  n  4 E(0  s)  4X*(8  s  +  t)  +  l  (l-cOct3  l  4  (6  xs) 4y*(t)  , 
i*0  A  j=o  i=0  e  * 


where  4e  (s)  and  4x*(s)  are  respectively  Laplace-Stieltjes 
transforms  of  e  defined  at  (9)  and  the  exponential  variable 
X*.  Setting  s  =  t  =  1  in  (15)  gives 


4X*  v*  (1»1)  =  E{exp(-X*)  exp(-X*  .)} 
*h'*n-k  n  n_i 


E<VW 


<1«> 


Then  using  the  fact  that  for  a  uniform  random  variable 
E(X)  ■  1/2  and  var(X)  =  1/12,  we  have  from  (15)  and  (16), 
after  simplification 


\  -  corr(X  ,X  k)  -  — ?  g  I  ( - ^ - r) 

k  n  n  k  2  +  eK  i-1  V  1  +  (l-oOe1/ 


k  *  1,2, •••  • 
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(17) 


Note  that  this  is  not  simply  a  geometrically  decaying  corre- 
lation  sequence,  as  for  the  NEAR(l)  process.  However,  for  the 
important  special  case  when  6=1  we  get 


and  thus  the  serial  correlations  are  the  kth  power  of 

p^,  which  takes  on  any  value  between  0  and  1.  Thus  we  have 
a  particularly  simple  uniform  Markovian  sequence. 

A  similar  analysis  given  in  Lawrance  and  Lewis  [5] 
shows  that 

E(Wl  -  n)  -  T  HI  +1  (1-1)6  >  (1  (l9» 

so  that  the  regression  is  not  linear. 

This  uniform  sequence  could  form  the  basis,  via  a 
probability  integral  transform,  of  many  other  sequences  with 
given  marginals.  However,  marginal  transformations  do  not 
preserve  correlation  structure,  as  shown  at  (17) ,  and  it  is 
therefore  useful  to  see  whether  sequences  with  marginals 
other  than  exponential  can  be  generated  from  (8)  for  special 
cases  with  a  suitable  choice  of  innovation  sequence  en<  This 
will  result  in  a  simple  process  with  autoregressive  Markovian 
structure . 
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4.  MARKOVIAN  SEQUENCES  WITH  SOME  OTHER  MARGINALS 

Although  an  exponential  distribution  is  a  common  assump¬ 
tion  for  positive  random  variables  met  with  in  problems  in 
operations  research/  it  is  too  narrow  an  assumption  to  encompass 
real  situations.  Therefore  parametric  distribution  models  are 
invoked  which  include  the  exponential  as  a  special  case  and 
which  allow  for  the  modelling  of  data  which  has  greater  or  lesser 
dispersion  than  exponentially  distributed  data.  Two  commonly 
used  models  are 

(i)  the  Gamma (k, A)  distribution  whose  probability  density 
function  is 


f/..,  _  AUxJ^e"** 

f(X)  *  - TO - 


k  >  0;  X  >  0;  x  >  0/ 


where  r(k)  is  the  complete  gamma  function,  and 
(ii)  the  (convex)  mixture  of  exponential  random  variables 


-A.x  -A,x 

f  (x)  =  ir^Aje  A  +  (l-ir^e  *  , 

0  <  <  X2;  x  >  0,  0  <  <  1  .  (21) 


The  Gamma  distribution  has  dispersion,  measured  by  the  coef¬ 
ficient  of  variation  C(X)  =  o(X/E(X)),  which  is  greater  than 
the  exponential  value  of  1  if  k  <  0  and  less  than  1  if  k  >  1. 
The  mixed  exponential  always  has  C(X)  2  1»  the  equality 
occurring  when  the  special  case  of  an  exponential  random  vari¬ 
able  with  parameters  A^  or  X2  holds. 
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4A.  The  Gamma  GAR(l)  process 

Direct  solution  of  equation  (6)  using  Laplace-Stieltjes 
transforms  gives  [1]  that,  in  the  stationary  sequence,  for  the 
X  to  be  Gamma (k, A)  we  must  have 


♦e  (s) 


E(e'se)  = 


p  +  (1-p) 


T+s 


(22) 


For  k  integer  this  has  an  explicit  inverse.  For  example, 

2 

for  k  -  2  the  innovation  e  is  zero  with  probability  p  , 

is  exponential (X)  with  probability  2p(l-p)  and  is  Gamma (2, A) 

2 

with  probability  (1-p)  .  It  is  easy  to  show  in  general  that  e 

V 

is  zero  with  probability  p  ,  so  that  the  "zero  defect"  is  not 
serious  for  large  k.  A  method  of  simulating  a  random  variable 
whose  Laplace-Stieltjes  transform  is  equation  (22)  was  derived 
by  Lawrance  [7],  using  the  fact  that  this  sequence  arises  in 
a  particular  type  of  shot  noise  process.  Thus  we  have  the 


Gamma  Innovation  Theorem 

Let  N  be  a  Poisson  random  variable  with  parameter 
0  =  -k  &n(p).  Let  U^,  Uj , . • .  ,  UN  be  uniformly  distributed 
over  (0,1)  and  independent.  Let  Y. YN  be  exponential (A) 
and  independent.  Then  c  can  be  simulated  using 

N  u 

e  *  I  Ypm  if  N  >  0  , 
m-1  m 

=  0  if  N  =  0  .  (23) 
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A  proof  is  not  given  here.  Note  that  e  is  zero  with  prob- 
ability  exp{-k  £n(p)}  =  p  .  Also  the  Poisson  number  N 
of  uniform  and  exponential  random  variables  which  must  be 
generated  for  each  e  has  expected  value  0  =  -k  £n(p). 

This  will  be  prohibitively  large,  and  the  simulation  will  be 
very  inefficient,  if  k  is  large  and/or  p  is  close  to  zero. 
Neither  of  these  cases  is  serious,  however.  If  k  is  large, 
say  greater  than  50,  the  sequence  is  almost  normal  and  the 
usual  normally  distributed,  AR(1)  linear  process  can  be  used. 
If  p  is  as  small  as  0.001  then  E(N)  is  still  only 
k  x  (6.9078)  which  is  still  reasonable.  However,  for  p 
this  small  the  sequence  is  approximately  i.i.d.  Gamma  and 
acceptance-re jection  techniques  for  simulating  Gamma  variables 
are  known. 

It  is  quite  simple  to  adapt  Algorithm  1  to  the  GAR(l) 
case.  It  would  pay  to  have  a  built-in  routine  for  generating 
the  Poisson  variable  which  will  bypass  further  calculations 
if  N  *  0.  In  other  words  routines  for  generating  Poisson 
variates  which  start  by  searching  at  the  median  of  a  table 
of  cumulative  Poisson  probabilities  will  be  inefficient. 

Unfortunately  the  NEAR(l)  process  does  not  appear  to 

extend  to  the  Gamma  case;  it  can  be  shown  explicitly  that 

there  is  no  innovation  eR  in  equations  (8)  and  (9)  which 

will  make  X_  have  a  Gamma  distribution  with  k  *  2. 
n 
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4B.  Mixed  Exponential  Markovian  Process 

Fortunately  first-order  autoregressive  Markovian 
processes  with  mixed  exponential  marginal  distributions  can 
be  obtained  from  equations  (8)  and  (9)  in  two  special  cases, 
and  these  sequences  should  be  widely  useful  in  modelling 
stochastic  systems. 

(i)  The  case  a  =  1;  MEAR(l) . 

In  (1)  it  is  shown  that  the  solution  to  the  Laplace 
transform  of  eR  for  the  linear  model  (6)  is  a  constant  p 
plus  a  (generally)  non-convex  mixture  of  three  exponential 
functions.  This  can  be  shown  to  be  a  proper  density  function 
if  p  £  but  it  can  also  be  shown  that  it  is  not  a 

density  function  for  all  p  less  than  one  and  greater  than 
or  equal  to  zero.  However,  Lawrance  [6]  showed  that  unless 
is  much  smaller  than  X2  (and  thus  the  Xn  are  very 
over-dispersed  relative  to  an  exponential  random  variable) 
a  solution  exists  for  for  all  p.  Thus  we  have  a  use¬ 

ful  process,  although  again  the  zero-defect  of  order  p  is  a 
problem. 

(ii)  The  case  8  ■  1;  MEAR(l) . 

When  8=1  in  equation  (8),  a  mixed  exponential  process 
TMEAR(l)  is  obtained  which  is  extremely  simple  to  simulate 
since  the  innovation  en  is  just  the  mixture  of  two  exponentials 
for  all  0  £  p  <  1.  Moreover,  the  process  has  no  zero-defect. 
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As  discussed  above,  the  sample  paths  will  tend  to  "run  up," 
but  this  is  no  great  problem  unless  p  is  fairly  large. 
Thus  we  have  the  following  Theorem  which  we  state  without 
proof: 


TMEAR ( 1 )  Theorem 

Let  the  first-order  autoregressive,  Markovian  sequence 

{X  >  be  defined  by 
n 


- 


en  +  VnVl  ' 


n  *  1,2,3,. .. 


where  P{V  =1}  =  1  -  P{v  =0}  =  a  for  0  <  a  <  1.  Then  the 
n  n  — 

sequence  {XnJ  is  stationary  and  has  a  (convex)  mixed  exponential 
marginal  distribution  with  probability  density  function 

-A.x  -X9x 

fx(x)  =  ir^X^e  +  (l-ir^e  z  , 

0  <  <  X2j  0  <  it  <  1;  x  >  0  ,  (24) 


if  en  is  i.i.d.  and  has  a  mixed  exponential  distribution  given 
by 


“Ylx  "y2x 

fe(x)  »  OjYje  +  n2r2e  » 

Yi  >  y2  >  0;  >.  0;  n2  =  l-r^  , 


(25) 


where 
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ni  *  a  ”  VYjJ/tY  1  -  Y2>  ; 

Yx  2  *  (b  +  (b2  -  4ac)1/2}/2a  ; 

a  =  (l-a)y1y2; 
y  =  ir2y2  +  =  E(X)  ; 

b  =  y2  +  yx  -  ay  ; 

8  =  y2  +  yx  -  y;  y2  *  1/yx;  w2  =  1^Y2 

and  Xq  is  independent  of  e^,  ...  and  has  probability 

density  function  (24) . 

Note  that  the  special  cases  where  tt^  =  0  or  ir2  =  1 
give  NEAR(l)  exponential  processes  with  parameters  X2  and 
X^  respectively.  Thus  they  should  be  handled  by  Algorithm  2 
since  they  will  cause  computational  problems.  The  case 
X^  *  X2  also  gives  a  NEAR(l)  process  and  is  excluded  for 
similar  reasons. 

Another  computational  problem  arises  from  the  fact 
that  the  probability  density  function  (25)  for  e  is  not  a 
convex  mixture  if,  as  is  possible,  is  greater  than  one. 

Of  course,  if  <  1  then  e  is  generated  as  an  exponen¬ 
tial  (y^)  random  variable  with  probability  and  as  an 

exponential (y2)  with  probability  (l-n^)  .  In  the  other  case 
we  use  the  following  theorem  (see  e.g.  (1)): 
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Simulation  of  Nonconvex  Mixed  Exponential 

Let  5  =  ^  ^1“Y2A1)  •  Then  with  probability  1-tS 
the  innovation  e  is  an  exponential (y^)  random  variable 
and  with  probability  5  the  innovation  e  is  the  sum  of 
an  exponential (y^)  random  variable  and  an  exponential (y2) 
random  variable. 

It  would  be  useful  to  have  a  mixed  exponential  solu¬ 
tion  for  the  sequence  (8)  for  all  8  for  convex  mixed 
exponential  marginal  distributions,  but  this  seems  difficult 
to  prove  because  of  the  algebra  involved. 


5.  GENERALIZATIONS 


In  all  of  the  processes  discussed  here  the  correla¬ 
tions  are  non-negative  and  geometrically  decreasing.  A  scheme 
for  obtaining  alternating  correlations  which  are  possibly 
negative  is  given  in  [1]  and  [5],  Another  problem  is  that 
different  types  of  dependence  and  higher-order  Markovian 
dependence  might  be  encountered  in  data.  Schemes  for  obtain¬ 
ing  mixed  autoregressive  moving  average  exponential  sequences 
where  the  autoregression  has  order  p  and  the  moving  average 
has  order  q  are  given  in  (4).  The  mixed  exponential  process 
TMEAR(l)  is  easily  extended  to  give  a  process  with  this  type 
of  extended  correlation  structure.  This  will  be  discussed 
elsewhere . 
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